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Abstract 

In this paper we propose a numerical method to solve the Cauchy problem based 
on the viscous shallow water equations in an horizontally moving domain. More 
precisely, we are interested in a flooding and drying model, used to modelize the 
overflow of a river or the intrusion of a tsunami on ground. We use a non conservative 
form of the two-dimensional shallow water equations, in eight velocity formulation 
and we build a numerical approximation, based on the Arbitrary Lagrangian Eule- 
rian formulation, in order to compute the solution in the moving domain. 
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1 Introduction 



The flooding and drying of a fluid on ground is a problem which has several 
applications in fluid mechanics such as coastal engineering, artificial lake fill- 
ing, river overflow or sea intrusion on ground due to a tsunami. These works 
involves the modelling of the physical process near the triple contact line be- 
tween the liquid (the water), the gas (the atmosphere) and the solid (the 
ground). The numerical simulation of this problem is very complex, due to 
the domain deformation and the triple contact line movement. 

Some numerical approaches to solve this problem have been previously ex- 
plored. For example, O. Bokhove |3] uses a conservative form of the shallow 
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water equations associated to the discontinuous galerkin discretization. D. 
Yuan et al. [16] propose to use a non conservative form of the shallow water 
equations in total flow rate formulation and uses a finite difference scheme for 
the numerical computation. 

In this paper, we use a two dimensional viscous shallow water model in eight 
velocity formulation. This model, called Stokes-like by P.L. Lions ([6], section 
8.3) is, in some sense, intermediate between the semi-stationary model and 
the full model of compressible isentropic Navier-Stokes model. This model 
can be obtained by integrating vertically the so-called primitives equations 
of the ocean [S] with some hypothesis on the viscous terms. We have chosen 
this model because we can prove the existence of a solution under certain 
considerations as the smoothness of the initial conditions and an acceptable 
hypothesis on a boundary operator [8]. 

After the presentation of the model we describe the numerical method based on 
the demonstration of the convergence proposed in [8] . Finally, we present some 
numerical examples in two idealized domains and we give some perspectives 
for other studies. 

1.1 Mathematical model 

We assume that we know a continuous function H from to M that represents 
the topography (the level of the ground with respect to a reference level). 

We note VLt (an open simply connected set of M?) the horizontal domain oc- 
cupied by the fluid at time t and rj{t, x, y) the elevation of the fluid compared 
to a horizontal zero level. We set /i(t, x, y) the eight of the water column, 
h{t^ a;, y) = H{x, y) + ri{t, x, y). 

The shallow water equations are based on a depth integration of an incom- 
pressible fluid conservation laws in a free surface-three dimensional domain. 
Governing equations for u and h can be obtained in the usual way (pj for 
example) and the two-dimensional system can be written as follows [5|7|^ : 

i'Au + Cd\u\u = fi, in y {t} X Qt (1) 

tG[0,T] 

in U {t}xnt (2) 

te[o,T] 

where g is the gravity constant, u the viscosity coefficient, a drag coefficient, 
/i represents external forcing (for example the wind stress). These equations 
are completed by initial conditions {u(t = 0) = Uq, rj{t = 0) = rjo) and 
boundary conditions. In the following, we assume that the boundary of the 
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real domain is non- vertical. Hence the variation of the surface elevation at the 
boundary imposes a modification of the horizontal domain occupied by the 
fluid. 



1.2 Boundary conditions 

1.2.1 In a three dimensional domain 

In many geophysical applications, the governing equations are solved in a 
domain [0,T] x Vi^ whereby VL^ = {{x,y,z) e R^/{x,y) G fit, —H{x,y) < 
z < 7]{t,x,y)} is a three dimensional moving domain. We denote by = 
(rii, n2, n^) the exterior unit normal to this domain, and we set Nh = {ni, 77-2) as 
its horizontal component. Moreover we denote V = (fi,f2,f3) as the velocity 
of the fluid and Vh = (f 1, V2) as its horizontal component. 

In the following, we assume that the ground of the sea does not have a vertical 
wall. Then we have ^ 0. The impermeability boundary condition on the 
ground for the three dimensional problem imposes V{x, y, —H{x, y)) ■ N = 0. 
Then, we have Vh ■ = —v^n^ on the ground and, since 7^ 0, we obtain 
V3 = -Vh ■ Vif. 




Fig. 1. Vertical section of the three dimensional domain 

At the surface, vertical velocity is equal to the lagrangian derivative of the 
function t], V3 = ^ + Vh ■ Vr]. 

We recall that we are interested in the effect of the overflow and we need to con- 
sider the problem in a moving domain that represents the domain actually oc- 
cupied by the fluid. Then we need to impose a condition on the moving bound- 
ary. At the triple contact line (between air, ground and water, see figured]), 
the horizontal velocity at the surface {vh{x, y, rj, t)) is equal to the velocity at 



3 



the bottom {vh{x,y, —H,t)). Similarly, we have V3{x,y,r],t) = vz{x,y, —H,t), 
coherent with the fact that the lagrangian derivative of h is equal to zero on 
this boundary. 

This property is not verified if, during its displacement, a part of the triple 
contact line touches a vertical wall. In this case, the boundary cannot be 
represented by the line h = 0. This question is very interesting but not treated 
in this paper. 



1.2.2 In a two dimensional domain 

The mean velocity u used in the equations [1] and [2] is related to the previous 
three dimensional velocity hj u = {ui, U2) = j!lfj{vi,V2)dz where h 0. The 
two dimensional domain fit is the projection of the three dimensional domain 
on the horizontal plane. 

For the depth integrated model, we denote by 70 the boundary of Qq (at 
initial time). Assuming that 70 is smooth enough, we define the deformed 
boundary as follows : 74 = {x' = x + d{x,t), x G 70}, where d represents the 
horizontal displacement d{x, t) = x(0, t,x)—x and where t, x) denotes the 
Lagrangian flow, i.e. the position at time t of the particle located at x at time 
s. With these notations, we have the natural condition dd[x{^,t,x),t)/dt = 
U{x, t) on 7(, where U{x, t) = u{x + d{x, t),t). 

The continuity equation ([2]) can be written as follows: 

h{x{xo,0,t),t) = h{xo,0)exp divv^ . (3) 

So, if at the initial time, h{xo,0) = then at all times /i(x(a;o, 0, t), t) = 0. 
Moreover, if h{xo,0) 7^ then /i(x(xo, 0, t), 0) 7^ at all times. We conclude 
that the boundary (and only the points of this boundary) at the initial time 
can define the boundary at all times. For example, if at initial time the domain 
is simply connected, the domain remains simply connected at all times t ^ 0. 

We characterize the boundary motion 74 by a condition on the normal com- 
ponent of the fluid stress tensor a on the boundary 

a{x + d{x,t),t) ■ n{x + d{x,t),t)\detJ\{x,t) = A (^—^f^ > (4) 

where A is a differential operator defined on 70 taking into account the bound- 
ary forcing on the fluid (see [2] for example) and J is the jacobian matrix asso- 
ciated to the transformation x h- )■ x(2;, 0, t). In [8], the authors show that if the 
operator A is a Laplace-Beltrami operator {j^^A{v)v = J^^ A^^'^(v)A^^'^{v) = 
||'y||^2(^Q)) and if the data are small enough, the problem has a solution. The 
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condition on A is necessary to obtain the smoothness of the boundary and 
allows the use of the classical Sobolev injections [HITS]. More precisely, in 
this paper, D{A) = if^(7o) with p = 2 and currently, we are not able to prove 
an existence result for p < 2. But numerically, it is possible to obtain physical 
acceptable results with less smoothness on this operator (for example in [7] 
we have obtained results by taking A = 0). 

It is not easy to characterize the operator A. It needs a good physical inter- 
pretation and give sufficient mathematical smoothness to prove the existence 
of a solution of our problem. Physically, this operator represents a friction 
condition in the triple contact line (between the fluid, the solid domain and 
the air). In the literature, a large number of physical approaches in the case 
of a wetting film of fiuid or for a droplet [Ij are found. Unfortunately, we 
cannot use directly these results. Firstly, we do not work at the microscopic 
scale (mainly governed by the Van Der Waals forces) and we need to take into 
account the effect of these phenomena at large scale. Secondly, the main effect 
studied for those phenomenon is in the vertical direction (for example for the 
capillary effect) but here, we use a depth integrated model. 



2 The numerical method 



2.1 The ALE equations 



We use a numerical method based on the weak formulation of the problem in 
the domain Qt depending on time. The ALE method that here we use, has 
been used to solve the Navier-Stokes equations in a moving domain [15]. 

We assume that at the initial time to, the fluid domain is covered by a regular 
mesh (for example, an usual finite element mesh). We also assume that on the 
boundary, each point of this mesh has the velocity of the fiuid. The velocity of 
the internal points can be chosen arbitrarily, the only condition is to conserve 
the smoothness of the mesh at all times. 

Then, at any time t G [0, T], the velocity c* of the moving mesh can be defined 
by solving the following problem: 



Ac* = in 

c* = u{x, t) on 7t 



(5) 
(6) 
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and we consider the following mapping 

{x,t) h-). c{x,t) = c*(x), 

where Q = Ute[o,T]{^t x {t})- The relation between the different domains fi^- 
is given by the mapping 

c{;ti,t2) : nt, ^Vlt, 

Xi^ X2 = C{xi,ti,t2) 

where C(.,ti,t2) is the characteristic curve from (xi,ti) to (x2,t2) correspond- 
ing to the velocity c* in the space-time domain Q. For each time r, we denote 
by Ut- and h^- the ALE velocity and the ALE thickness: 

Ur{x,t) = u{C{x,T,t),t), hr{x,t) = h{C{x,T,t),t). (7) 

With these notations, we can write the ALE formulation of the problem 
Ou 

+ {Ur — Cr)'VUr — flAUr + Crf Ur\Ur\ + k^. = 0{t — t) , (8) 

dh 

-—^ + (Ur-Cr)VK + hrdiYUr=0(t-T). (9) 

ot 



In what follows, we will use a first order time discretization to solve this 
problem (see [10] for more details about second order schemes). 



2.2 Time discretization 



Let M ^ 1, we note At = as the time step. For 1 ^ n ^ M, 1 ^ m ^ M, 
we set fi" = Qt^ with boundary and 

glix)=gtAx,t'')iorxen^, 
g-{x)=g:^ix)=g{x,t'^)ioTxEn\ 

In order to ensure the positivity of ht^ , the continuity equation is renormalised 
as follows 

d 

— log hfn + (ut^z — Qn) V log ht^ + div Ut^ = 0. (10) 
dt 

With this renormalisation, we do not have conservation of the mass, but when 
we follow the evolution of the mass during a simulation (e.g. in the first test 
presented in the following section), only small variations of the mass quantity 
are observed. 
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We denote by f/" (respectively and C") the approximation of the exact 
solution -u^n (respectively /i^n and Qn). 

Then, we note f/"(x) = f/"(X"(x,r)) and R^{x) = i7'^(X"(x, T)) with 
X"'{x, ■) the characteristic curve solution of: 

'^(x,t) = (f/"-C")(X"(x,t)) 

_X"(x,r+^) = x. 

In our study, we approximate the foot of the characteristic by X"'(a;, t") ~ 
x—{U"'—C'^){x)At. We set U'^~^^ (resp. -ff^"^^) the approximation of Utn{x, 
(resp. /;,fn(x, t""*"^). Approximating the Lagrangian derivative by a first order 
Euler scheme, and using a linearised drag operator, we obtain the following: 

z/Af/;^+i) = [7", infi" (11) 
in Q"", (12) 

on 7", (13) 
on 7". (14) 



log + At div f/^+^ = log if" 



divf/, 



n+1 



A 



dt 



curl f/' 



■n+1 



h: 



n+1 



Implicitly, taking into account the equation [31 we have H^~^^ = on 7". This 
is an approximation at the first order in time of the initial shallow water 
problem. 

Remark 2.1 

This problem is solved by a fixed point technic. A first approximation of\l^ is 
computed assuming U^~^^ = [/" and the approximation of H^~^^ is then used in 
[77] to obtain a first approximation ofU^~^^. We repeat this operation in order 
to obtain convergence of the system. Prove of convergence can be fount in JB]. 

We cannot easily increase the order of this scheme because the ALE formula- 
tion is only of the first order. 

The previous problem is solved on by using the spatial discretization pro- 
posed in the following section. After, we need to solve the problem given by 
equations [5] and [6] in order to compute the mesh velocity. Each point pk with 
coordinates Xk of the mesh is then moved with the first order approximation 
x^+i =a;^ + AtC"(a;^,r). 



7 



2.3 Spatial discretization 



The spatial discretization of the previous problem is based on the Finite El- 
ement Galerkin method. In [5], an approach based on the Galerkin method 
with a special basis is proposed, but this approach cannot be applied easily if 

We note {0^} the set of finite element functions of H^{Vr). The weak formu- 
lation of our approximated problem is 



gAt H:^\x)diYM^) + CdAt U:+\x)\U:ix)\M^) 



+ / a{X + diX,t),t) ■ n{X + diX,t),t)^{Mx)) 
J "ft 

= I (15) 
\ogHl+\x) + Atdivf/^+i(x) = log#"(x) in fi" (16) 

where the operator 7 is the trace operator from to 7". 
Then the equations (jlj) and f[T^ give 

/ U:+\x)M^) + uAt f VU:^\x)VM^) 

-gAt [ H:+\x)divM^)+cJ U:^\x)\U:{x)\M^) 

f/"(x)0fc(x) (17) 



where g{x) = g{x{x,0,t)). 

We then use a first order discretization of the partial time derivative: 



A'/' (^) (7(0^(^0))) ^70 - ^ A^" (f/r^(x)) A^" {liMxo))) djo 
1 1 A'/' (tj^ix)) A'/' (7(^(^0))) ^70. (18) 



The first part of the right hand of this equation needs to be included on the 
left hand of the global problem and the second part on the right hand. 
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The global weak formulation is: 

JCl" JQ" JQn 

Ur\x)\U:(x)\Mx) + /^^ (c/r^(x)) A'/' (7(0.(xo))) rf7o 

= J^^ U-{x)Ux) + A''' {i{k{xo))) d^o, (19) 

log H^+^ (x) + At div C/;^+^ (x) = log (x) . (20) 

Remark 2.2 On jt, U{l,t) = J2kC'k4>k{f',t) where I is a point of jt (md the 
sum is computed on all finite element functions. Then, on the initial boundary 
7o, we have the decomposition lJ{lo,t) — J2k'^k<Pk{lo,'t) "^iih the same Uk- 

Remark 2.3 Formally, we can use two kinds of boundary conditions: 

(1) A condition of normal displacement of the boundary. We write u-n = ^ 
where d is the normal displacement. With this condition, the well posed 
weak formulation associated to the diffusion operator is /QdivMdiv0 + 
/QCurlucurl0 — /^divM0-n + f^cmlucj)- a{n) where a{ui,U2) = {—U2,ui) 
and curl u = ^ — ^ . The boundary condition is then imposed on h — 

ay ax ^ 

div u and curl u since we do not want to impose a condition on u ■ n or 
u ■ a{n) . 

(2) A condition on the displacement of the boundary u = ^ where d is a 
vector. With this condition, the well posed weak formulation associated to 
the diffusion operator is VuVcf) — 'S/u(f)-n and the boundary condition 
is then imposed on Vu (more precisely, we need to take into account the 
trace tensor including condition on h in the boundary, but we assume 
h = on this boundary). 

Remark 2.4 We can note that the friction term Cciu\u\ is necessary to sta- 
bilize the flow. Indeed, if we use the following decomposition 

{L\n'') f = WHliST) ® Curli/^(Q") © WH{n''), (21) 

where "H is the intersection of and the space of harmonic functions, we can 
see that the functions ofVHiVt^) are not controlled by the laplacian operator. 
We present in the following section a case with = where we do not have 
damping of the flow. 

2.4 First numerical test 

The previous method is tested in order to simulate the behavior of a fluid in 
a simplified domain. We assume that the domain is axisymmetric using the 
function H{r) : r i->- ar^ + 1. At the initial time t — 0, the fluid at rest touches 
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the wall {h{a,t = 0) = 0) (see figure [2]). Physically, our initial domain Qq is a 




Fig. 2. Simplified domain - (a) at rest, (b) at an arbitrarily time. 



disc with a radius of 130 meters. This domain is meshed with triangles (420 
triangles for the mesh Ml, 470 for the mesh M2 and 1002 for the mesh M3). 
The height of the column of water at the centre is 1 meter (since r = at the 
centre of our domain). Fluid viscosity is O.Olm^.s"^ and gravity coefficient g 
is assumed equal to lm.s~^ in order to amplify the elevation. 

We apply a forcing at the surface of this fluid. This forcing is usually taken 
into account by applying the continuity of the horizontal stress tensor on the 
surface. So, for the three dimensional model. 



where Wa is the wind velocity (at ten meters over the flow for the ocean for 
example) and C is a "drag coefficient" assumed to be constant. Using the 
vertical integration of the vertical operator -§^i^z^§^, and assuming Uh{Ti) — u, 
we obtain a forcing condition on the mean velocity / oc ^"jj^"^ . 

We use here 



Hence, we observe an oscillation of the free surface of the water plan. At each 
oscillation, a part of the water moves on ground and a part of the ground 
is uncovered by this water (see figure [3]) • In the last part of the remark 12. 4^ 
we indicate that with C^^ = a part of the flow component is not diffusive. 
More precisely, our solution is only a gradient. If we do not take into account 
boundary friction effects or bottom friction effects, a part of the fluid is not 
diffusive. To observe this effect, we plot (on flgure H]) the level of a bound- 
ary point according to time. We can observe an accumulation of energy on 




(22) 
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Fig. 3. Oscillation of the water 




Fig. 4. Oscillations of the water without drag term 

the fluid due to the numerical approximation, and this accumulation is not 
compensated by diffusive term. 

In figure |5l we plot the variation of the same point taking into account the 
drag coefficient Cd\u\u. With this term, all the modes of the fiuid are diffusive 
and, if we stop the forcing, global energy of the fluid vanishes. 

In flgure [6], we plot the time evolution of the kinetic energy of the flow for 
the different meshes (Ml, M2 and M3). We only have a little difference for 
all these meshes but a very expensive cost for the mesh M3. After the forcing 
phase (20 seconds), we can observe the transformation of the kinetic energy in 
potential energy and conversely. When the kinetic energy vanishes, potential 
energy is maximal. If the drag coefficient is assumed to be equal to zero, the 
amplitude of the oscillation does not decrease. 

Finally, figure [7] represents the evolution of the mass of fiuid compared to the 
initial mass. Even if we do not have a rigorous conservation of the mass, due 
to the renormalisation of the mass equation, the variation of the mass is very 
little. 
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Oscillation OT a boundary point 
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Fig. 5. Oscillations of the water with a drag term 



Evolution of the kinetic energy 
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Fig. 6. Evolution of the solution for Ml, M2 and M3 meshes 



Estimation of the mass conservation 
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Fig. 7. Evolution of the mass in comparison with initial mass 
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2.5 Second numerical test 



In this second numerical experiment, we use a more complex domain to test 
our numerical method. We assume that the domain is axisymmetric using the 
function H{r) : r i— i- ar^ + + or + d. b,c and d are chosen in order to at the 
distance 1/a of the centre of the fluid domain, the fluid touches the wall (see 
figure [8]). We use the same external forcing as for the first experiment fl22|) . 




Fig. 8. Mesh of the ground and initial position of the water 

But, due to the specific form of the topography, a part of the fiuid goes to the 
external crown. 

We present some results of the simulation in figure M Recalling that the shal- 
low water equations are based in the continuum mechanic, it cannot cut this 
domain with a finite energy. Hence, even if the thickness of the layer of wa- 
ter is very low, we conserve a thin film of water between all the parts of the 
domain occupied by the fiuid. We recall that the computation is only made 




Fig. 9. Evolution of the fluid domain 
in the moving two-dimensional domain Qf Mathematical study proves that 
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the main difficulty is to conserve the smoothness of the boundary. This result 
can be observed in this simulation because, at the final time, we have contact 
between two parts of the boundary (see last part of the figure [H]). Since we do 
not use a "good" boundary operator, we do not have sufficient smoothness on 
the boundary. We observe then a change of connectedness and a hole appears 
in the fiuid domain. 



3 Conclusion 

The presented work suggests that capillarity effect needs to be incorporated 
into the evolution equations, and more precisely in the triple contact line. 
Theoretical results presented in [8] give sufficient smoothness for the capillarity 
term used in equation HI 

The main difficulty is to describe the capillarity effects on the depth inte- 
grated model and in a dynamical system. A possible approach is to test some 
boundary operator and to compare numerical and experimental result. This 
approach will be proposed in future studies. 
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